Skip to contents
library(findSVI)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)
library(leaflet)
library(knitr)
library(tidyr)
library(glue)
library(htmltools)
library(sf) #not suggests
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

Get census data with geometry

First, we are using get_census_data() to obtain ZCTA-level data with simple feature geometry for PA for 2020 (Census API required).

pa_zcta_2020_geo_data <- get_census_data(year = 2020, state = "PA", geography = "zcta", geometry = TRUE)

Here, we are showing the first 10 rows of the data. With the geometry = TRUE, we’ll get a tibble with an additional column containing simple feature geometry (MULTIPOLYGON).

#> Simple feature collection with 10 features and 134 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -80.48268 ymin: 40.10397 xmax: -79.73383 ymax: 40.83308
#> Geodetic CRS:  NAD83
#> # A tibble: 10 × 135
#>    GEOID NAME        B06009_002E B06009_002M B09001_001E B09001_001M B11012_010E
#>    <chr> <chr>             <dbl>       <dbl>       <dbl>       <dbl>       <dbl>
#>  1 15001 ZCTA5 15001        1655         417        5450         582         724
#>  2 15003 ZCTA5 15003         441         125        2352         245         250
#>  3 15004 ZCTA5 15004          62          43          63          70          21
#>  4 15005 ZCTA5 15005         293          92        1380         200          65
#>  5 15006 ZCTA5 15006          37          58          13          28           0
#>  6 15007 ZCTA5 15007           0          11          22          35           0
#>  7 15009 ZCTA5 15009         404         121        3253         384         206
#>  8 15010 ZCTA5 15010        1237         242        5636         444         719
#>  9 15012 ZCTA5 15012         638         180        2750         302         336
#> 10 15014 ZCTA5 15014         136          75         527          96         105
#> # ℹ 128 more variables: B11012_010M <dbl>, B11012_015E <dbl>,
#> #   B11012_015M <dbl>, B16005_007E <dbl>, B16005_007M <dbl>, B16005_008E <dbl>,
#> #   B16005_008M <dbl>, B16005_012E <dbl>, B16005_012M <dbl>, B16005_013E <dbl>,
#> #   B16005_013M <dbl>, B16005_017E <dbl>, B16005_017M <dbl>, B16005_018E <dbl>,
#> #   B16005_018M <dbl>, B16005_022E <dbl>, B16005_022M <dbl>, B16005_023E <dbl>,
#> #   B16005_023M <dbl>, B16005_029E <dbl>, B16005_029M <dbl>, B16005_030E <dbl>,
#> #   B16005_030M <dbl>, B16005_034E <dbl>, B16005_034M <dbl>, …

Get SVI with geometry

After getting the data ready, we can supply this tibble with simple feature geometry to get_svi().

pa_zcta_2020_geo_svi <- get_svi(2020, data = pa_zcta_2020_geo_data)

get_svi() will return the full SVI table with every SVI variables, intermediate percentile ranks, theme-specific SVIs and overall SVI (consistent with CDC/ATSDR SVI database, without MOE).

pa_zcta_2020_geo_svi %>% head(10) %>% kable()
GEOID NAME geometry E_TOTPOP E_HU E_HH E_POV150 E_UNEMP E_HBURD E_NOHSDP E_UNINSUR E_AGE65 E_AGE17 E_DISABL E_SNGPNT E_LIMENG E_MINRTY E_MUNIT E_MOBILE E_CROWD E_NOVEH E_GROUPQ EP_POV150 EP_UNEMP EP_HBURD EP_NOHSDP EP_UNINSUR EP_AGE65 EP_AGE17 EP_DISABL EP_SNGPNT EP_LIMENG EP_MINRTY EP_MUNIT EP_MOBILE EP_CROWD EP_NOVEH EP_GROUPQ EPL_POV150 EPL_UNEMP EPL_HBURD EPL_NOHSDP EPL_UNINSUR EPL_AGE65 EPL_AGE17 EPL_DISABL EPL_SNGPNT EPL_LIMENG EPL_MINRTY EPL_MUNIT EPL_MOBILE EPL_CROWD EPL_NOVEH EPL_GROUPQ SPL_theme1 SPL_theme2 SPL_theme3 SPL_theme4 RPL_theme1 RPL_theme2 RPL_theme3 RPL_theme4 SPL_themes RPL_themes
15001 ZCTA5 15001 MULTIPOLYGON (((-80.43758 4… 31129 16070 14093 5567 787 3037 1655 1319 6395 5450 4871 815 142 5764 955 537 94 1189 619 18.2 4.8 21.5 7.0 4.3 20.5 17.5 15.9 5.8 0.5 18.5 5.9 3.3 0.7 8.4 2.0 0.5313 0.5595 0.5174 0.4001 0.4804 0.5516 0.3258 0.5877 0.6929 0.6555 0.8173 0.7801 0.4422 0.5305 0.7071 0.8031 2.4887 2.8135 0.8173 3.2630 0.5028 0.7418 0.8173 0.8224 9.3825 0.8059
15003 ZCTA5 15003 MULTIPOLYGON (((-80.2366 40… 11212 6084 5104 2656 212 1234 441 534 2090 2352 1656 356 86 1771 196 38 8 587 26 23.7 3.5 24.2 5.5 4.8 18.6 21.0 14.8 7.0 0.8 15.8 3.2 0.6 0.2 11.5 0.2 0.7039 0.3761 0.6615 0.2798 0.5417 0.4240 0.5840 0.4991 0.7989 0.7469 0.7860 0.6724 0.2872 0.3886 0.8234 0.5585 2.5630 3.0529 0.7860 2.7301 0.5380 0.8473 0.7860 0.6515 9.1320 0.7707
15004 ZCTA5 15004 MULTIPOLYGON (((-80.38652 4… 380 207 138 75 0 6 62 38 42 63 61 21 0 30 0 25 0 0 0 19.7 0.0 4.3 22.5 10.0 11.1 16.6 16.1 15.2 0.0 7.9 0.0 12.1 0.0 0.0 0.0 0.5837 0.0000 0.0490 0.9393 0.8597 0.0959 0.2690 0.5997 0.9641 0.0000 0.6118 0.0000 0.7499 0.0000 0.0000 0.0000 2.4317 1.9287 0.6118 0.7499 0.4745 0.2520 0.6118 0.1266 5.7221 0.2287
15005 ZCTA5 15005 MULTIPOLYGON (((-80.24356 4… 9191 4275 3948 916 123 609 293 250 2534 1380 1288 92 19 459 134 23 0 135 213 10.2 2.5 15.4 4.1 2.8 27.6 15.0 14.3 2.3 0.2 5.0 3.1 0.5 0.0 3.4 2.3 0.2295 0.2593 0.2234 0.1896 0.2890 0.8678 0.1986 0.4594 0.2963 0.5312 0.4801 0.6684 0.2775 0.0000 0.3613 0.8212 1.1908 2.3533 0.4801 2.1284 0.0999 0.4597 0.4801 0.4631 6.1526 0.2707
15006 ZCTA5 15006 MULTIPOLYGON (((-79.88531 4… 292 131 131 0 0 0 37 0 58 13 0 0 0 0 0 0 0 0 0 0.0 0.0 0.0 13.3 0.0 19.9 4.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0000 0.0000 0.0000 0.7911 0.0000 0.5091 0.0545 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.7911 0.5636 0.0000 0.0000 0.0392 0.0148 0.0000 0.0000 1.3547 0.0136
15007 ZCTA5 15007 MULTIPOLYGON (((-79.93657 4… 629 167 167 387 0 51 0 19 414 22 295 0 0 0 0 0 0 0 0 61.5 0.0 30.5 0.0 3.0 65.8 3.5 46.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9920 0.0000 0.8627 0.0000 0.3186 0.9881 0.0499 0.9903 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2.1733 2.0283 0.0000 0.0000 0.3785 0.2985 0.0000 0.0000 4.2016 0.1022
15009 ZCTA5 15009 MULTIPOLYGON (((-80.45619 4… 15114 7057 6509 1904 323 1424 404 315 3718 3253 2196 271 57 869 528 18 23 335 304 12.8 4.5 21.9 3.6 2.1 24.6 21.5 14.8 4.2 0.4 5.7 7.5 0.3 0.4 5.1 2.0 0.3109 0.5259 0.5436 0.1635 0.2186 0.7832 0.6317 0.4991 0.5322 0.6226 0.5176 0.8199 0.2479 0.4479 0.5105 0.8031 1.7625 3.0688 0.5176 2.8293 0.2236 0.8593 0.5176 0.6850 8.1782 0.6118
15010 ZCTA5 15010 MULTIPOLYGON (((-80.48268 4… 27205 12699 11363 4761 811 2322 1237 923 5815 5636 4584 820 42 3819 811 251 135 1243 1268 18.3 6.1 20.4 6.6 3.4 21.4 20.7 17.1 7.2 0.2 14.0 6.4 2.0 1.2 10.9 4.7 0.5359 0.7077 0.4536 0.3678 0.3759 0.6101 0.5647 0.6547 0.8114 0.5312 0.7611 0.7954 0.3806 0.6632 0.8046 0.9030 2.4409 3.1721 0.7611 3.5468 0.4784 0.8939 0.7611 0.9222 9.9209 0.8610
15012 ZCTA5 15012 MULTIPOLYGON (((-79.90417 4… 15243 7709 6883 2639 376 1186 638 470 3595 2750 2629 348 86 1048 353 445 22 481 144 17.3 5.0 17.2 5.7 3.1 23.6 18.0 17.2 5.1 0.6 6.9 4.6 5.8 0.3 7.0 0.9 0.4972 0.5806 0.2917 0.2991 0.3311 0.7412 0.3570 0.6581 0.6279 0.6924 0.5681 0.7356 0.5436 0.4137 0.6370 0.7003 1.9997 3.0766 0.5681 3.0302 0.3099 0.8627 0.5681 0.7491 8.6746 0.7015
15014 ZCTA5 15014 MULTIPOLYGON (((-79.74922 4… 3055 1671 1450 761 41 334 136 175 464 527 379 162 0 264 48 0 7 136 91 26.0 2.5 23.0 6.2 5.9 15.2 17.3 12.8 11.2 0.0 8.6 2.9 0.0 0.5 9.4 3.0 0.7625 0.2593 0.6057 0.3331 0.6599 0.2287 0.3121 0.3481 0.9242 0.0000 0.6328 0.6575 0.0000 0.4752 0.7533 0.8587 2.6205 1.8131 0.6328 2.7447 0.5573 0.1947 0.6328 0.6578 7.8111 0.5409

For visualization purposes, we’ll simplify the table and keep only the GEOID, geometry and SVIs.

svi <- pa_zcta_2020_geo_svi %>% 
  select(GEOID, geometry, contains("RPL_theme"))
svi %>% head(10) %>% kable()
GEOID geometry RPL_theme1 RPL_theme2 RPL_theme3 RPL_theme4 RPL_themes
15001 MULTIPOLYGON (((-80.43758 4… 0.5028 0.7418 0.8173 0.8224 0.8059
15003 MULTIPOLYGON (((-80.2366 40… 0.5380 0.8473 0.7860 0.6515 0.7707
15004 MULTIPOLYGON (((-80.38652 4… 0.4745 0.2520 0.6118 0.1266 0.2287
15005 MULTIPOLYGON (((-80.24356 4… 0.0999 0.4597 0.4801 0.4631 0.2707
15006 MULTIPOLYGON (((-79.88531 4… 0.0392 0.0148 0.0000 0.0000 0.0136
15007 MULTIPOLYGON (((-79.93657 4… 0.3785 0.2985 0.0000 0.0000 0.1022
15009 MULTIPOLYGON (((-80.45619 4… 0.2236 0.8593 0.5176 0.6850 0.6118
15010 MULTIPOLYGON (((-80.48268 4… 0.4784 0.8939 0.7611 0.9222 0.8610
15012 MULTIPOLYGON (((-79.90417 4… 0.3099 0.8627 0.5681 0.7491 0.7015
15014 MULTIPOLYGON (((-79.74922 4… 0.5573 0.1947 0.6328 0.6578 0.5409

Interactive maps for overall SVI

With the simple feature geometry, we can visualize SVI patterns and perform spatial SVI analysis with any mapping tool. Here, we’re using a powerful package leaflet for interactive maps.

First we’ll examine the missing value.

missing <- svi %>% filter(is.na(RPL_theme1)) 
missing %>% kable()
GEOID geometry RPL_theme1 RPL_theme2 RPL_theme3 RPL_theme4 RPL_themes
15260 MULTIPOLYGON (((-79.95539 4… NA NA NA NA NA
15616 MULTIPOLYGON (((-79.563 40…. NA NA NA NA NA
15691 MULTIPOLYGON (((-79.68701 4… NA NA NA NA NA
16312 MULTIPOLYGON (((-79.30513 4… NA NA NA NA NA
17077 MULTIPOLYGON (((-76.53944 4… NA NA NA NA NA
17120 MULTIPOLYGON (((-76.88651 4… NA NA NA NA NA
17822 MULTIPOLYGON (((-76.60949 4… NA NA NA NA NA
18225 MULTIPOLYGON (((-75.97559 4… NA NA NA NA NA
18430 MULTIPOLYGON (((-75.49468 4… NA NA NA NA NA
18936 MULTIPOLYGON (((-75.23686 4… NA NA NA NA NA
19109 MULTIPOLYGON (((-75.16434 3… NA NA NA NA NA
19112 MULTIPOLYGON (((-75.19692 3… NA NA NA NA NA
19442 MULTIPOLYGON (((-75.58142 4… NA NA NA NA NA

Looks like there are a few ZCTAs with missing values in SVIs. (Upon checking the full SVI table, we can see that most of the individual SVI variables are 0 for those ZCTAs.) So we’ll remove those ZCTAs for better visualization.

To set up the interactive map for overall SVI, leaflet() does most of the heavy lifting. Here, we’ll just add some customized color palette and labels for the appearance.

svi_clean <- svi %>% drop_na()

#set color palette
pal <- colorNumeric(
  palette = c("orange","navy"),
  domain = svi_clean$RPL_themes
)

#set label
zcta_label <- glue("<h3 style='margin: 0px'>{svi_clean$GEOID}</h3>
                    overall SVI: {svi_clean$RPL_themes}") %>%
  map(~HTML(.x))

leaflet(svi_clean) %>% 
  addProviderTiles(providers$CartoDB.Voyager) %>% 
  addPolygons(color = "white", 
              weight = 0.5,
              smoothFactor = 0.5,
              opacity = 1,
              fillColor = ~pal(RPL_themes),
              fillOpacity = 0.8,
              highlightOptions = highlightOptions(
                weight = 5,
                color = "white",
                fillOpacity = 0.8,
                bringToFront = TRUE),
              label = zcta_label,
              labelOptions = labelOptions(
                style = list(
                  "font-family" = "Fira Sans, sans-serif",
                  "font-size" = "1.2em"
                ))
              ) %>% 
  addLegend("bottomleft",
            pal = pal,
            values = ~RPL_themes,
            title = "Overall SVI in all ZCTAs in PA (2020)",
            #labFormat = labelFormat(prefix = "$"),
            opacity = 1)
#> Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
#> Need '+proj=longlat +datum=WGS84'

From the interactive map, we can visualize easily how SVI varies in different regions in PA and zoom in to examine specific ZCTAs of interest, making it a helpful approach to explore new ideas, patterns and analyses.